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Abstract. We present a numerical study of the mechanical response of a 2D Lennard- Jones 
amorphous solid under steady quasistatic and athermal shear. We focus here on the evolution 
of local stress components. While the local stress is usually taken as an order parameter in 
the description of the rheological behaviour of complex fluids, and for plasticity in glasses, we 
show here that the knowledge of local stresses is not sufficient for a complete description of 
the plastic behaviour of our system. The distribution of local stresses can be approximately 
described as resulting from the sum of localized quadrupolar events with an exponential dis- 
tribution of amplitudes. However, we show that the position of the center of the quadrupoles 
is not related to any special evolution of the local stress, but must be described by another 
variable. 
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1 Introduction 

While the mechanical plastic response of crystalline 
solids is now well described in terms of elementary 
defects called dislocations, the plastic response of 
amorphous solids is not so well understood, and 
has been the object of numerous studies in the past 
two decades. Most theoretical works [I|[2||3l[4l[5lf6l[7l 
11111111111^^ . on the forced 

dynamics of glassy materials assume that the ir- 
reversible evolution of glasses occurs through lo- 
calized, irreversible events. Such events have been 
well identified in simulations of 2D models [I] and 
are referred to as shear transformation zones (STZ) 
or quadrupolar events 18, 19 , 20 , 21 and arc some- 
what similar to Tl events in the rheology of 2D 
foams [22] . 

These localized and dissipative events should 
correspond to the "soft mode" of the material [TBI 
H9 l [23 j [24] ) . which is the first mode to become unsta- 
ble, when the system reaches a limit of mechanical 
instability. However, the local and intermittent na- 
ture of the plastic events complicates the identifica- 



tion of the corresponding order and control param- 
eters. Experimental studies on disordered materi- 
als, far below the glass transition temperature [22j 
[Z3fcMnKMMM2IX!H\, support the existence of 
a collective behaviour of localized rearrangements 
leading to a strongly heterogeneous mechanical re- 
sponse as shown in the existence of shear bands 
in the macroscopic plasticity of such systems. It 
has also been shown experimentally [32 j [33 | f34] as 
well as numerically [il[l^pTp5|[36l[37] while study- 
ing the stress-strain response of amorphous ma- 
terials that, in the so called plastic regime, these 
events lead to an intermittent release of the macro- 
scopic stress in the material. The understanding 
of the role played by local stresses in the macro- 
scopic mechanical behaviour of such systems is thus 
very important. Since these systems are known to 
have a wide distribution of local stresses, meso- 
scopic models [71IT0 1 ITT 1 I38 ] propose precisely to re- 
late the existence of shear bands to the local rear- 
rangements due to local thresholds and activated 
dynamics in the local stress itself. In these mod- 
els, the collective effects would thus be due to long- 
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range elastic couplings in solid systems, and the ir- 
reversible behaviour would be associated with local 
stress thresholds, and avalanche dynamics. 

In the present paper, we provide an extensive 
study of the local stress dynamics in the quasi-static 
plastic response of a model glass at very low tem- 
perature, and in the plastic flow regime. The sys- 
tem we have studied is a 2D Lcnnard- Jones glass, 
already studied previously for understanding the 
heterogeneous mechanical response in the elastic 
regime, far below plastic flow [55 | l40 | l4i p 42 | l2T]. In 
the first part of this paper, we discuss the evolu- 
tion of the local stresses (local pressure, deviatoric 
stress) in the plastic flow regime, within the frame- 
work of dynamical heterogeneities [33] ■ We show 
the existence of a finite size of cooperativity ap- 
pearing in the plastic flow regime, even in the limit 
of zero temperature and quasi-static shear. In the 
second part of this paper, we provide a detailed 
analysis of the full statistical properties of the stress 
components (pressure, shear stress). We show that 
the temporal correlations shown in the local parti- 
cle rearrangements for small imposed strain [211119) 
are absent in the average evolution of local stress 
components. We then show in the third part of this 
paper, that the statistical evolution of stress com- 
ponents can be described as resulting simply from a 
sum of localized quadrupolar rearrangements. This 
description allows us to propose a simple equation 
of evolution for the local stress depending unfortu- 
nately on a local yet unknown criterion for instabil- 
ity. In the last part, we show that this local criterion 
for instability is not related to any stress threshold 
but could be related to a different variable, such as 
the local stress increments. 



2 Dynamical heterogeneity 

2.1 Two points correlation function 

The system we study numerically has already been 
described extensively in previous articles [4T? l l4T | l4"2" l 
[21"]). It is a 2D glass obtained by quenching very 
quickly a liquid sample obtained by molecular dy- 
namics simulations, and made of polydisperse par- 
ticles interacting through a Lennard- Jones poten- 
tial. In the following, all quantities will be expressed 
in Lennard- Jones units (ulj), that is an energy 
clj = 1(« O.leV), a distance olj = 1(~ 2A) and 
mass m = 1(~ 10g. mo? -1 ). The number of parti- 
cles per unit surface is p — 0.925, corresponding to 
very small pressure P w 0.2 at rest. The system is 
made of at least N = 10000 particles, correspond- 
ing to a square with a linear size of 104 interatomic 



distances. Larger systems have also been consid- 
ered, especially in this paper a system of N = 23125 
particles, corresponding t a rectangular box of size 
50 x 500. These large system sizes are necessary in 
order to avoid strong finite size effects. The quench- 
ing starts from the liquid state, at a temperature 
T = 2. Then the system is relaxed using MD to a 
given temperature of 1 — 0.5 — 0.1 — 0.05 — 0.01 — 
0.005 - 0.001, where it is aged during 1000 unit 
times. The final step is a relaxation at zero tem- 
perature, using the Conjugate Gradient Method in 
order to get an instantaneous minimization of the 
potential energy. The system being at zero tempera- 
ture, it is sheared quasi-statically, by imposing suc- 
cessive steps of global shear strain 5e = 5.10 -5 on 
the walls. After each step, the system is relaxed into 
its nearest equilibrium position, thus ensuring that 
the applied mechanical deformation is quasi-static. 
The local stress components are computed on each 
particle i by using the usual Irving-Kirkwood for- 
mula ([33]) 

<Ta0,i = —1/Vi.Ejtijda-dp.rij 

Vi being the volume of the Voronoi cell, Uj the 
amplitude of the interacting force between parti- 
cles i and j, the distance between particles and 

d the unit vector of the bond. The macroscopic 
stresses are obtained by averaging the local stress 
components over the whole system, or alternatively 
by measuring the force per unit length acting on 
the walls along appropriate directions. In Fig.l, we 
show the (averaged) macroscopic shear stress as 
a function of the total strain applied to the sys- 
tem. The global shear stress has an intermittent 
behaviour, with an alternance of small and large 
jumps, giving rise to negative slopes that are the 
signature of dissipative events. Note that the dis- 
tribution (figure 2) of these incremental stresses is 
peaked around the value given by the elastic re- 
sponse with an elastic shear modulus \i « 11.7 ([40)) 
and shows around this value softer, as well as more 
rigid, steps. 

In order to characterize the evolution of the stress 
in the flowing regime, we first compute the two 
points correlation function of the local stress com- 
ponents acting on each particle inside the sample. 
The correlation function is "Lagrangian" in the sense 
that the stress is followed on a given particle as the 
system is sheared. In the plastic flow regime, we 
expect the system reaches a stationary steady state 
with timc-translational invariant (TTI) observables 
and correlations. Note that in quasistatic simula- 
tion, time does not appear explicitly and must be 
interpreted in the following as the number of steps n 
of elementary shear strain 8e applied to the system. 
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The total imposed strain is thus e = n.Se. In order 
to study the stress fluctuations correlation function, 
we introduce the notation p(r , n) = cr(V , n) — o~(n) 
where A denotes a spatial average and "F* is the po- 
sition of the particle; a is a component of the stress 
tensor, namely a = Odev where <7dev = &i — 02, <J\ 
and (T2 being the eigenvalues of the local 2x2 stress 
tensor; or a = a p where a p = — (a± + 0"2)/2 is the 
local pressure; or a = is the local shear stress. 
The two-time autocorrelation function of the stress 
fluctuations writes 



successive exponential decays (see the fit by an ex- 
ponential decay with characteristic strain e w 10% 
at very large strains in Fig. 3) . In any case, the pres- 
ence of dissipative events leads to the relaxation of 
local stress at large deformations, above a few per- 
cents of strain. 



2.2 Four points correlation function or 
cooperativity number 



C(An, n) — p(r ,n+ An)p( r , n), (1) 



and 



C{An) = (C(An,n)} 



(p( r , to + An)p( r 
C(An,n), 



(2) 



where the last equality assumes TTI. 

Numerically we have computed the stress com 
ponents (as mentioned before from the Irving-Kirkwoogerature; for larger applied strains, in the plastic 



In the recent years, the quest for a dynamical coop- 
erativity length associated with the slowing down of 
the dynamics of supercooled liquids [43 l l47 ( [48] has 
led to the development of new statistical tools to 
characterize such "dynamical heterogeneities". In 
our previous paper |21j we showed that the motion 
of the particules in our system under shear is highly 
non trivial. It shows a background of heterogeneous 
motion even in the elastic regime at very low tem- 



definition, denoted by o~i(n)) on each particle of the 
system. The discrete version of the above equation 



writes 



C(An) = 



1 

TV 



Piin + An)pi{n)). 



(3) 



where Pi(n) = <7j(n) — a(n). We have computed 
this quantity in the plastic flow regime (i.e. for 
external shear deformation e greater than 1.5%), 
where the average total stress saturates and fluctu- 
ates around a given average value cry. We checked 
that the correlations are invariant under time trans- 
lation and in Fig. 3 we plot the autocorrelation 
function for the shear component of the stress, on 
a lin-log scale and normalized to 1 for An = 0. The 
behaviour can be decomposed into two successive 
relaxations. It starts with a first exponential decay 
at small deformation with a characteristic strain 
of 3.1%, this decrease to a plateau is analogous to 
the P — relaxation in supercooled liquids, and is 
followed by an apparent stretched exponential de- 
cay. A similar Kohlrausch- Williams- Watts [45]l46] 
is usually observed in supercooled liquids and age- 
ing glassy systems, with a correlation that evolves 
as C(t) = Aex-p(-{t/r)P) with f3 < 1. This be- 
haviour is approximately recovered in our systems 
at large deformations (see inset of Fig. 3), with the 
exponent j3 = 0.7. For small strains, the relaxation 
is strictly exponential. Note that the large strains 
stretched exponential decay could also be seen as 
the best fit for an intermediate regime between two 



regime, it shows zones of very high mobility located 
in the vicinity of elementary shear bands and of lo- 
calized quadrupolar events, and similarly zones of 
low mobility far from these irreversible events (the 
local displacement field can vary by many orders 
of magnitudes). This type of behaviour in the me- 
chanical response of glass, even at zero tempera- 
ture, can be seen as some kind of dynamical hetero- 
geneity (while the "dynamics" here is overdamped) 
and it is of great interest to quantify a cooperativ- 
ity degree or a cooperative length scale. To do so, 
Bcrthier ct al. [43] have proposed to look at the so 
called Xi f° ur point correlation function. An other 
possible observable is the cooperativity number in- 
troduced by Doliwa and Heuer [47] , that estimates 
the spatial fluctuations of the two-point correlation 
function introduced previously. Following reference 
l47l the definition we use here is 



N rj 



Var[SX t ] ({^^-(Ei^)} 2 ) 



SVar[Xi 



(4) 



whereA 4 (t',t) = {a^t+t')-^^ + t'))(a % {t)-a t {t)) 
is a dynamical quantity associated with every parti- 
cle i, and computed on various components a of the 
stress tensor. In Eq.(4), one shows easily that in the 
case of totally uncorrelated variables Xi ((XiXj) = 



0), one gets N c , 



1, but in the opposite if the cor- 



relation is complete (Xi = Xj) one gets N coop = N. 
Finally in the case of L independent groups each 
comprising M identical variables and of zero mean 
(TV = LM) one gets 
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3 Statistical analysis of the local shear 
stress 



N C oop 1 



= l 



N M-1 



i 



= 1 3 = 1 



The cooperativity number is obtained by av- 
eraging over the origins t and remains a function 
of time t' (N coop (t')). Note that this function does 
not contain any useful information on the heteroge- 
neous motion in the linear elastic regime, since in 
this case, the amplitude of the stress on each parti- 
cle increases linearly with time. The cooperativity 
number is thus constant in this regime. A maxi- 
mum can be reached only in the non-linear regime, 
like in the plastic flow regime here, or when ther- 
mal activation starts playing a role. In Fig. 4 we 
have plotted this function computed on all parti- 
cles located at a distance of a few (typically 10) 
particle diameters from the boundaries, and for dif- 
ferent components of the stress tensor. All the stress 
components have approximately the same contribu- 
tion to the cooperativity number, with a maximum 
of a few tens at an imposed strain of a few per- 
cents (s=s 4%). As described before, the maximum is 
a measure of the maximum number of particles on 
which the local stress evolves in a correlated man- 
ner. In our case, it corresponds to an ensemble of 
6 x 6 to 7 x 7 particles where the stress components 
evolve in a cooperative way. The strain at which 
this maximum is reached is about 4 pcrcents, that 
is beyond the plastic threshold, and close to the 
characteristic strain obtained from the two points 
correlation function C. This picture corresponds to 
the description of localized plastic events, that we 
have made already in a preceding paper |21| . In the 
clastic regime (as in elastic branches present in the 
plastic flow plateau [21]), the cooperative evolution 
of stress can be more important, but does not con- 
tribute significantly to the cooperativity number. 

Fig. 4 shows that there are strong boundary ef- 
fects measurable on this quantity. Indeed, the co- 
operativity number is much higher in the vicinity 
of the walls (« 100) corresponding to a layer of 
one particle diameter along the entire system size 
(L x = 104). Finally, note that this maximum coop- 
erativity number is finite in the quasi-static limit 
studied here, suggesting a saturation at very low 
shear rate. 



In the previous section we have studied the average 
evolution of the stress components by looking at 
averaged quantities like the number of cooperativ- 
ity, or the two-point correlation function. We will 
now examine in details the full histogram of stress 
changes. The results are shown here for the shear 
stress component, but we have checked that all the 
following results are also valid for other components 
(pressure and deviatoric stress) . In figure 5 we show 
the typical evolution of the shear stress component 
on a given particle up to a total strain of 259o 
The variations of stresses are larger than expected 
if they would be Brownian. In order to analyze this 
evolution, we have plotted in figure 6 the distribu- 
tion P(Aa xy , An) of the shear stress increments 

Aa xy (An, n) = ( o xy (n + An) — a xy (n + An)) 
- (vxy(n) -W X y(n)) (6) 

for various numbers An of incremental shear steps 
(An = 1,2,4,8,16...), averaged over the origins n 
and over the entire system. The distributions P(Aa xy , 
are all symmetric, with zero mean. For An — > oo 
one recovers a gaussian distribution, which is con- 
sistent with the central limit theorem. However, at 
small imposed strains, these distributions are not 
gaussian as would be the case for a Brownian evo- 
lution. A finer analysis of the distribution of the 
elementary increments P(Aa xy , An = 1) (see fig 
6b and c) shows three zones. At small incremental 
stress jumps Aa xy , we can see a plateau of approx- 
imately constant probability, whose width evolves 
inversely proportional to the volume V of the sam- 
ple, followed by an apparent power-law decay, in a 
zone of approximately three decades for 100/1/ < 
Aa xy < 10 where P(Aa xy , An = 1) oc 1/At£+\ 
and concluded by an exponential cut-off indepen- 
dent on the system size (characteristic shear stress 
Aa xy w 1.4). This upper cut-off allows for a finite 
variance of the local stress evolution. In the absence 
of any temporal correlations, the entire process can 
thus be described by the central limit theorem, with 
a scale invariant distribution on a finite interval. In 
this intermediate stress range, the probability den- 
sity function (hereafter, pdf) P(Aa xy , An) is well 
reproduced by the scale invariant relation 



1 In a previous paper we have analyzed in a similar 
way the motion of a particle in the system under shear 



(see Ref. [21]) 
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Aa. r 



P(Aa xy ,An) = An~ H f(^f) 



d(a(i,e) - cr(e)) 
de 



(9) 



Ai 
Aa 



P(Aa xy ,Ae) = Ae-«f(^f) 



with 



f(u) oc 



-a-l 



, for u <C 1 
, for 1 <C u <C cste.An~ H 



(7) 



(8) 



In the intermediate stress range, the process can 
thus be considered as self-similar. Fig 7(a) illus- 
trates this scaling with a good superposition of the 
distributions for a — 0.7 and H ~ — (only for not 
to large An since for very large An the upper cut-off 
discussed before contributes significantly to the re- 
sulting distribution). The exponent a describes the 
algebraic (slow) asymptotical decay of the distri- 
bution of the incremental jumps P(Aa = s, An = 
1) oc s~ a ~ 1 (as shown in fig 6(a)). The exponent 
H is related to the evolution of the stress jumps 
as a function of the applied shear strain An. It 
characterizes the Z\n-dependence of the crossover 
between a regime of approximately uniform proba- 
bility (for \Aa\ <C An H ) and the power law regime 
(for \Aa\ ^> An H ). The coefficient H also accounts 
for possible temporal statistical correlation between 
jumps. As stated in Taqqu et al. [49] the only non 
degenerate a — stable self similar processes with sta- 
tionary increments that verify H = — and where 
< a < 1 are the a — stable Levy motions. As 
described above, the evolution of the shear stress is 
thus of Levy flight type, but only in the interme- 
diate stress (and applied strain) range. The Levy 
flight evolution implies by definition, first that there 
is no temporal correlation between local stress jumps 
during the shear of the sample, second that the 
variance of elementary changes is infinite (as long 
as the cxponentiel cut-off is neg lectedfl These re- 
sults can be compared with the study of the pdf 
P(Ay, An = 1) of the positional jumps in the trans- 
verse direction, that showed in contrary that these 
jumps were correlated in time for small imposed 
deformation [21] . 

According to the above, a plausible equation 
that would describe the evolution within a Levy 
flight process for the stress component a, averaged 
over the whole system, is 



2 In the case of Brownian motion the properties stated 
above on the pdf and on the process are verified (self 
similarity, a — stable process, stationary increments) but 
with a finite variance (a > 2) allowing to a unique value 
H = 1/2. 



where n(i,e) is a stochastic process whose spatial 
average 77(e) is the process that is entirely charac- 
terized above (ie by the distribution of the elemen- 
tary increment and by the absence of correlation 
between successive increments), and e is the exter- 
nal imposed shear strain. For small imposed shear 
strain, it can be mentioned that such an equation 
with a noise corresponding to a Levy motion can- 
not be reinterpreted in terms of the usual Fokkcr- 
Planck equation as the moments of the shear stress 
a are non vanishing for all orders. However, for large 
strain intervals Ae, the stochastic process becomes 
gaussian, due to the existence of the upper cut-off 
in the distribution of rj. 

We will now propose a microscopic origin for 
the elementary process 77, by comparing it with a 
spatial redistribution of quadrupolar stress. 

But before, it can be interesting (as in Ref. [HO 
[20] ) to identify the contribution to the distribution 
P(A<t, An) of events with a release of the macro- 
scopic stress (plastic events with Aa macro < 0). 
This is shown in figure 7(b). We see here that for 
large imposed strains (An > 64 it means Ae > 
0.32% - smaller than the observed plastic thresh- 
old) the contribution to the distribution is mainly 
due to plastic events meaning that the weight of 
elastic events is negligible in this range of applied 
strain. This corresponds also to the strain interval 
An at which the distribution becomes gaussian with 
a good approximation (non gaussian parameter less 
than 1), and where the departure from Levy mo- 
tion is of course noticeable. For smaller strains, the 
contribution of steps with an increase of the macro- 
scopic stress (Admacro > 0) influences the distribu- 
tion mainly for small values of the local stresses. 
These positive steps are for example responsible for 
the different power-law decay (oc l/Aa^y in place 
of oc 1/Aa^y) in the beginning of the first cross- 
over in the distribution (at the end of the initial 
plateau). We are unable to explain the precise value 
of these two exponents. However, we will now pro- 
pose a model that is able to explain the generic fea- 
tures of this distribution of shear stress increments. 



4 Minimal model: sum of quadrupoles 

As shown in [21] , the plastic flow regime consists in 
a succession of elastic branches and plastic events. 
We will now focus on the contribution of plastic 
events to the evolution equation of the stress com- 
ponents in the system. We identify here a plastic 
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event by the fact that the associated stress is re- 
laxed macroscopically (corresponding to an evident 
dissipation of energy shown by a negative slope in 
the stress-strain response of the system, see Fig.l). 
Some of them correspond to an isolated quadrupo- 
lar redistribution of local stresses around a well 
identified center inside the sample; other plastic 
events consists in an alignment of rotational struc- 
tures in the particle displacements giving rise to ele- 
mentary shear bands going through the system [21] . 
The two kinds of events coexist in the plastic flow 
regime, and the transition from one kind of plastic 
event to another is still matter of debate [TM3"3l 
[50ll51j. However, the isolated quadrupolar event is 
more frequent than the elementary shear band. In 
the linear initial regime for example, only the iso- 
lated quadrupolar events take place, while there are 
periods with only quadrupolar events in the plastic 
flow regime as well. 

In the following, we take the quadrupolar iso- 
lated event as the elementary building block for 
explaining the plastic deformation of the material, 
and we propose a simple model that describes the 
plastic deformation (as shown in the redistribution 
of stresses) as a sum of uncorrelated quadrupoles of 
random amplitude A. 

We first identify the distribution P(A) in our 
system. As shown by Picard et al. [52], a quadrupo- 
lar event involves in 2D a long-range redistribution 
of stresses, due to a local pure shear. The corre- 
sponding stress change is of the form. 



Aa xy (r, 9) 



A% , if r < r < r ma 
A , if < r < r 



(10) 



where r and r max are respectively the typical size 
of the quadrupole and the size of the system, A is 
the amplitude of the quadrupoles, and Aa xy de- 
notes the incremental shear stress. We have ne- 
glected in this expression the quadrupolar angular 
dependence of the stress field and only considered 
its p. spatial decay 0. The distribution of P(A) 
is measured in our data (see figure 8). It corre- 
sponds to the incremental stress at the center of 
the quadrupole, that is at the place where the dis- 
placement is maximum. The best fit is exponential, 
with P(A) = 9^j- ex P(~~ ~t~) an d &a = 2 appears to 
be the characteristic amplitude of the quadrupolar 
event. It is independent on the system size, unlike 
the distribution of macroscopic stress release that 
is oc 1/V (figure 2), and is symmetric. 



The angular dependence will mainly contribute by 
a scaling factor, and affect the weight of small Aa in a 
logarithmic way in the final result. See M. Tsamados et 
al. preprint (2007) for a detailed calculation. 



Using the exponential distribution P(A) shown 
above, we can now reproduce the distribution P(Aa = 
s, An = 1) of the incremental shear stress Aa xy 
averaged over the whole system. According to the 
previous analysis (see part III) , the distribution for 
An > 1 should then be reproduced by successive 
convolutions, assuming the absence of temporal cor- 
relations. Neglecting the angular dependence of Aa xy 
as in Eq. I10| we show that, for a given amplitude 
A, there is a simple bijective relation between the 
shear stress Aa xy and the radial coordinate r. We 
can thus write PA(Aa xy )dAa xy = P(r)dr, with 
P(r)dr = 2-Kr/Vdr. Using Eq[TDl we get 



P A (Aa xy ) 



irAr 2 
VAa 2 



with 



\A\- 



< \Aa xy \ < \A\ 



(11) 



(12) 



In the following, we will consider only the case A > 
0, the opposite case giving the symmetric distribu- 
tion corresponding to Aa xy < 0. Assuming that the 
successive events contribute independently to the 
total distribution P(Aa xy ) (an assumption justified 
by the fact that this distribution is averaged over 
the time origins) , the average probability of having 
an incremental stress Aa xy inside the system is thus 
obtained by summing over all quadrupolar events 



P(Aa xy ) = c 



P A (Aa xy )P(A)dA (13) 



with an additional normalisation factor c due to the 
limited range of allowed amplitudes. This gives 

7TT 2 

P ( Aa xy ) = c— ^ 2 {{a A + Aa xy )exp(~Aa xy /a A ) 

xy 



■'-)exp{~A<j xy .r 2 max /{a A .rl))) (14) 



x v r 2 
'o 

We show here that the upper exponential cut-off 
is proportional to exp(— Ao xy joA) and that the 
low Aa behavioi0 is dominated by finite size scal- 
ing proportional to exp(— Aa xy .r 2 lax / {(ja-Tq))- We 
have plotted this function in figure 9, and com- 
pared it with the result obtained numerically. We 
find a good agreement with the numerical result, 
for ga — 2, r max = L = 100 and r = 4. The pa- 
rameter ro is the unique free parameter of the fit 
since r max scales like the system size and the fitted 



4 with a logarithmic divergence if the angular depen- 
dence is taken into account see M. Tsamados et al. 
preprint (2007). 
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value for a a corresponds to the value obtained in 
the measured distribution of P(A) (see figure 8), 
thus confirming our simple model. The small value 
obtained for r shows that the plastic quadrupolar 
events are localized. 

This simple model gives us a physical intuition 
on the origin of the distribution for the increments 
rj observed in the averaged evolution equation [9l 
Therefore we can now make more precise the evo- 
lution equation at a local level as 



d(<r(?,e)-<r(e)) 
de 



d 3 i 



A{e)rl 



II r -r (e)|| : 



(15) 

where tq is the position where the quadrupolar 
event has occurred, || r — tq\\ > ro , and p(rQ,e) 
is the number of plastic event per unit strain and 
unit volume, taking place at tq when the external 
strain is e. 

As it appears clearly here, this equation is not 
complete since rg(e) is unknown. In many models 
[TTlllOp the occurrence of a quadrupolar event at 
position Fq is triggered by a local stress threshold. 
In the following section, we discuss the existence of 
such a criterion. 



5 Looking for a local threshold 

The center of the quadrupolar plastic rearrange- 
ments is identified as the place where the particle 
displacement is the largest [21] . In order to relate 
the position of the center of the quadrupoles to local 
stress criteria, we have compared the distribution of 
the stress components obtained for the whole sys- 
tem, and the value of the same stress component 
on the particle in the center of the quadrupole, one 
step before the quadrupole takes place. 

The result for the deviatoric stress adev = G \ ~ 
02 is shown in figure 10(a). The distribution of 
the deviatoric stress in the whole sample is sta- 
tionary: it is the same in the beginning (for ex. 
one step before the first plastic event) and in the 
end (e.g. one step before the last measured plastic 
event) of the plastic flow. It is remarkable to sec in 
this figure, that the distribution of the deviatoric 
stress restricted to the center of the quadrupoles 
(one step before each plastic event) computed for 
all the plastic events is very close to the previously 
discussed stationary distribution obtained over the 
whole sample, and not only in the center of the 
quadrupoles. The position of the center of the quadrupohmponents evolve in a heterogeneous way during 
appears to be unrelated with any threshold value in quasi-static plastic deformation. It is related to a 
the local stress components. The saturation of the finite cooperativity number identifying fluctuating 



macroscopic stress, and the associated well defined 
yield stress ay in these systems can thus be related 
only to the alternance of increase and decrease of 
stress and to its intermittent behavior, rather than 
to any identified local stress threshold. The exis- 
tence of an apparent macroscopic yield stress is 
not related to a local yield stress. Note that cri- 
terion based on Tresca or Mohr-Coulomb criterion, 
involving the most probable relation between devi- 
atoric stress and pressure have also been tested in 
our numerical systems. It shows a general tendency 
for the pressure to be affinely related to the devi- 
atoric stress in the plastic flow regime. However, 
this tendency is shared by all the particles in the 
system, and not only by the particles at the cen- 
ter of the quadrupoles. It thus appears, like in the 
case of the deviatoric stresses, that a global Mohr- 
Coulomb criterion, or equivalently -for 2D systems- 
a Tresca criterion, is valid in the plastic flow regime, 
but is does not provide a locally selective criterion 
for plasticity. 

Figure 10(b) shows the same distributions as for 
the deviatoric stress, but for the incremental devia- 
toric stresses: it means (7<2 e „(n+ 1) — Odevin)- In this 
case, quenched stresses are not taken into account, 
and a large incremental stress is generally the signa- 
ture of a large local deformation. The distribution 
of incremental stresses is also stationary. The dis- 
tribution restricted to the center of the quadrupoles 
is displaced to larger increases in the deviatoric 
stresses. This point has also been mentioned by 
Robbins et al. [12] and underlines the role played 
by incremental stresses in the dissipative dynamics 
of the systems, in comparison with total stresses 
that are more or less unchanged. It suggests that 
the incremental stress Aa can be a more important 
parameter for plastic purposes than the stress a it- 
self. However, its distribution is very broad, and 
a criterion based only on the incremental stresses 
is not very selective: it does not even exclude the 
possibility for the center of the quadrupole to take 
place where the increase of the deviatoric stress is 
not maximum but minimum... 

The question of the existence of a local crite- 
rion for plasticity appears thus clearly beyond our 
detailed numerical study of local stress components, 
and deserves a further studies. 



6 Conclusion 



We have shown in this paper that the local stress 
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zones whose maximum size is approximately 6x6 
particle diameters in the plastic flow regime. The 
strain at which these zones form (« 4%) is related 
to the typical strain at which the memory of the lo- 
cal stress seems to be lost in the 2-point correlation 
function. 

It is thus natural to propose a picture of the 
heterogeneous stress evolution, as a sum of local- 
ized plastic events. We have shown in this paper 
that the averaged evolution of the full distribution 
of local stress components can be described by a 
simple model. Within this model, the evolution of 
local stress components results from the uncorre- 
lated sum of localized quadrupolar rearrangements. 
The model allows to recover the 3 regimes shown 
in the distribution: its 1/V dependence at small 
stresses, in agreement with a local description of 
the plastic instability, the cross-over to an apparent 
power-low decay cx 1 /Aa a+1 with a close to 1 in the 
intermediate range, and for larger stresses a size- 
independent upper cut-off directly related in our 
model to the size-independent distribution of stress 
rearrangements in the center of the quadrupoles. 

However, a complete description of local stress 
rearrangements is still lacking: for example, the ob- 
served exponential nature of the distribution of the 
amplitudes of stress rearrangements in the center of 
the elementary quadrupoles - a crucial ingredient of 
our model - is still not explained, and the position 
of the center of the quadrupoles should be related 
to a criterion of instability not identified in this pa- 
per. This last question especially deserves further 
studies, since is appears clearly here that the local 
stresses are not determinant. 

Any attempt to describe local plasticity in terms 
of local stress rearrangements should thus take into 
account these numerical results before proposing a 
complete description of the heterogeneous mechani- 
cal response in terms of a single evolution equation. 
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Fig. 2. Histogram of the slope Aa/Ae of the macroscopic shear stress. Inset: Size dependence of the amplitude 
of negative stress changes Aa xy . The amplitude of stress releases is proportional to the volume of the sample. 
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Fig. 3. Autocorrelation function of the shear stress computed on each particle and averaged over the whole 
sample, as a function of time (applied shear strain). The nature of the fits is indicated in the legend box. Inset: 
Logarithm of the Autocorrelation function of the shear stress, in log-log scale, in order to determine the exponent 
of the corresponding stretched exponential 
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Fig. 4. (a) Cooperativity number as a function of time (applied shear strain) for various stress components: 
local shear stress a xy , local pressure P, eigenvalues of the stress tensor <n and CT2, and the deviatoric stress 
Odev — <f \ — °"2 (from bottom to top). The curves have been obtained by averaging on 20 000 time origins to. (b) 
Cooperativity number of the local shear stress at the borders, in comparison with the cooperativity number in 
the center of the sample 
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Fig. 5. Evolution of the shear stress on a given particle in the center of the sample, as a function of the applied 
shear strain. For each value, the average over the sample has been substracted. 
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Fig. 6. (a) distribution P(Aa xy , An) of the variation of local shear stress Ao xy for different strain intervals An. 
The distribution has been obtained by averaging over the whole sample, and for different time origins, (b) lin-log 
scale for P(Aa xy , An = 1) and different system sizes N = 10000(104 x 104) and N = 23125(50 x 500). It shows a 
size independent exponential upper cut-off with a characteristic Aa ~ 1.4. (c) log-log scale for P(Aa xy , An = 1) 
and different system sizes N = 10000(104 x 104) and N = 23125(50 x 500). It shows the scaling of the first 
cross-over oc 1/V where V is the volume of the system. The stress variation is mainly due to local variation that 
averages over the whole system (see text). 
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Fig. 7. (a) scaling of P(Aa xy , An) oc <f>{Aa xy / An 1 ^) with z = a and a = 0.7. The scaling is very good in the 
intermediate range of the distribution, and for small An < 64 where the exponential cut-off does not play a role. 
Inset: the same with a = 1 = z. This exponent corresponds also to the one discussed in the text, (b) separate 
contribution of steps associated with a release of stress at a macroscopic level, to the total distribution of local 
stresses, for time intervals An = 1 and An = 64 
. As shown here, the distribution of Aa for An > 64 is dominated by the contribution of stress changes due to 

plastic events. 
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Fig. 8. Distribution of shear stress releases in the center of the quadrupoles, during a plastic event. The fit is 
exponential with a characteristic Aa = 2. 
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Fig. 10. (a) Distribution of the deviatoric stress in a configuration of 10 000 particles for 2 different strains in the 
beginning and in the end of the plastic flow plateau (triangles). Comparison with the distribution of the deviatoric 
stress in the center of the future quadrupole just before the quadrupolar event takes place (squares), computed 
on the 2500 plastic events of the plastic flow regime. The 3 curves are very similar, (b) Same as in (a), but for the 
incremental deviatoric stress. In this case, the distribution obtained on the center of the quadrupoles indicates a 
more pronounced correspondance with the highest incremental stresses inside the sample. 



